
import random
import math

'''
    test maximum entropy principle
'''


def averager(r):
    return (r + 2*r**2 + 3*r**3)/(1+r+r**2+r**3)
    

def r_given_mu(mu):
    rl=0
    rr=10000
    avl=averager(rl)-mu
    avr=averager(rr)-mu
    assert avl*avr<0, 'something wrong with '+str(mu)
    rm=0.5*(rl+rr)
    avv=averager(rm)-mu

    while math.fabs(avv)>1e-7:
        if avv*avl>0:
            rl=rm
        else:
            rr=rm
        rm=0.5*(rl+rr)
        avv=averager(rm)-mu
        print rm, avv
    return rm



samples=[]
mu=1.43


parameter_r=r_given_mu(mu)
C=1./(1+parameter_r+parameter_r**2+parameter_r**3)
for i in range(4):
    print i, C*parameter_r**i
        
        
    
# nucleotides in genome
N=100


for i in range(100000):
    x=[]
    for i in range(N):
        x.append(random.randrange(4))
    
    sumx= sum(x)
    #print sumx
    if sumx==mu*N:
        samples.append(x[:])
        #print len(samples), '<< new!'
        

distr={}
for s in samples:
    for v in s:
        distr[v]=distr.get(v,0)+1

all_counts=sum(distr.values())
for gc, v in distr.iteritems():
    print gc, float(v)/all_counts

